#ifndef _CTRACKMATRIX_H
#define _CTRACKMATRIX_H
#include "CVector3d.h"
#include <string>
#include "public.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
using namespace std;
class CBeam;
class CTrackMatrix
{
	public:
		CTrackMatrix();
		~CTrackMatrix();
		CTrackMatrix(double,double ,int dh=1,int _flagE=1);
		CLHEP::HepMatrix &Drift(double,double _dEnergy=0);
		CLHEP::HepMatrix &Quad(double,double,double _dEnergy=0);
		CLHEP::HepMatrix &RFGap(double,double _dPhis=0,double _dEnergy=0);
		            void DtlCell(double,double,double, double,double,double,double,double,double _dEnergy=0);
								void DtlTank_File(string );
								void InitLattice();
								void Equip_Match();
			vector<double> Equip_DesignQuad(int _iperiod=0);
						     int Match_Conver(vector<vector<double> >&);
								 //vd2 is mismatch factor in x,y,z;返回值是三个方向最大的失配参数
							double Match_MMF(vector<double>& ,vector<double> &,int mt=12);
							  void Lattice_MMF(string );
							  void inSig_modify(double );
							   int Match_Follow(const vector<vector<double > >& vvd,int iStart=0,int mt=0);
							double Match(int nit=100,int iStart=0);//匹配成功返回0，
							double Match8(int nit,int iStart,int nc=4);
							double Match8_1(int nit,int iStart,int nc=4);
							   int nls_(int *nm, double *dd, int *mm,int iStart);
								void setDelt(double _dt=0.0056){Delt=_dt; };
								void test();
								//format: SL,SR,E0TL,Phi,L1,S1,L2,S2,EnergyIn
		            void BeamVector (vector<double>& );
		            void TwissVector(vector<double>& );
		CLHEP::HepMatrix TransTwiss(CLHEP::HepMatrix &);
		            void TransSig(CLHEP::HepMatrix &);
		CLHEP::HepVector TwissToEnve(CVector3d,CVector3d,CVector3d);
                void setFreq(double);
                void setCurrent(double);
                void setEnergyIn(double);
								void getTwiss();
								void getPhaseAdvance();
								void getSigMatrix();
                void outTwiss();
								void setPeriod(int i=1){Period=i;};
								void getStruct();
								void setFlagGapEmit(int i=0){flagGapEmit=i;};
			vector<double> CalSig(vector<CLHEP::HepSymMatrix>&,int);
			vector<double> CalSig3(vector<CLHEP::HepSymMatrix>&,CLHEP::HepMatrix &,int);
			//vector<double> CalSig(vector<CLHEP::HepSymMatrix>&,int,CLHEP::HepMatrix &,int _j=0);
			vector<double> CalSig2(const CLHEP::HepMatrix &hd);
                //Equipartitioning
                void Equip();
                void Equip_Tradtion();
								void Equip_Cell(double ,double,int _xflag=0,int _yflag=1);
								//void Equip_Cell(double ,double,double _zflag=0,int _yflag=1);
                //static double wrap(const std::vector<double> &x, std::vector<double> &grad, void *data);
  private:
    CLHEP::HepMatrix SC_linear(double);
	private:
    vector<CLHEP::HepVector>     vhdTwiss;
    vector<CLHEP::HepMatrix>     vhmatrix;
    vector<CLHEP::HepMatrix>     vhmatrix_Total; // store the matrix multiply till the current cell
    vector<CLHEP::HepMatrix>     vhmatrix_Period; // store the matrix multiply till the current cell
		vector<vector<double> >      vvdEmit;
		vector<vector<double> >      vvdBeta;
//		vector<vector<double> >      vvdBeamSize;
	public:
		CLHEP::HepMatrix    hepmatrix;
		CLHEP::HepVector    hepvector;
		CLHEP::HepVector    hepvectorTwiss;
		CLHEP::HepVector    beami;
		CLHEP::HepVector    beamf;
		CLHEP::HepVector    beamo;
		CLHEP::HepMatrix    hepmatrixTwiss;
		CLHEP::HepSymMatrix hepmatrixSig;
    double              dCurrent;
    vector<CLHEP::HepSymMatrix>  vhdSig;
		vector<vector<double> >      vvdStruct;
    vector<vector<double> >      vvdPh_Av;
    vector<vector<double> >      vvdPh_Av_Equip; //the lenth is the total number of periods
    vector<vector<double> >      vvdFF; //match
		vector<vector<double> >      vvdph_ad_0;//equip in 0 current
		vector<vector<double> >      vvdStruct_Equip_0;//equip in 0 current
	private:
		double dLQuad1,dSQuad1,dLeft,dRight,dLQuad2,dSQuad2;
		double dEnergyIn,dE0TL,dPhis,dFreq,dEnergyf;
		int    flagE,dh;  //flagE determines whether to count the energy gain in RFGap
		double dk0x,dk0y,dk0z,dkx,dky,dkz;
		double dEmitx,dEmity,dEmitz;
		int    Period,iCell;
		int    flagGapEmit;
		double  Delt;
//    double *dtest1;
};
#endif
